In [ ]:
import sys
globber_path = '/Users/adrian/projects/globber/'
if globber_path not in sys.path:
sys.path.append(globber_path)
from collections import OrderedDict
import astropy.units as u
import astropy.coordinates as coord
import numpy as np
from numpy.polynomial.polynomial import polyval
import matplotlib.pyplot as pl
import matplotlib as mpl
pl.style.use('apw-notebook')
%matplotlib inline
import h5py
from scipy.ndimage import gaussian_filter
from scipy.stats import scoreatpercentile
from scipy import interpolate
from globber.ngc5897 import cluster_c, r_c, r_t, fiducial_DM
In [ ]:
XCov_filename = "/Users/adrian/projects/globber/data/ngc5897/XCov_test.h5"
In [ ]:
with h5py.File(XCov_filename, "r") as f:
ra = f['search']['ra'][:]
dec = f['search']['dec'][:]
allX = f['search']['X'][:]
cluX = f['cluster']['X'][:]
nonX = f['control']['X'][:]
isoX = f['isochrone']['X'][:]
# print(list(f.keys()))
non_ll = f['log_likelihood']['control'][:]
iso_ll = f['log_likelihood/isochrone/15.60'][:]
In [ ]:
# all_data = np.array([tuple(x) for x in np.vstack([ra,dec,non_ll,clu_ll,xd_ll]+[il for il in iso_lls.values()]).T],
# dtype=[('ra','f4'),('dec','f4'),('non_ll','f4'),('clu_ll','f4'),('xd_ll','f4')] + \
# [("iso_"+name,'f4') for name in iso_lls.keys()])
all_data = np.array([tuple(x) for x in np.vstack([ra,dec,non_ll,iso_ll]).T],
dtype=[('ra','f4'),('dec','f4'),('non_ll','f4'),('iso_ll','f4')])
pre_filter_ix = (allX[:,0] > 18.3) & (allX[:,0] < 20.) & np.isfinite(all_data['non_ll'])
# pre_filter_ix = np.ones_like(allX[:,0]).astype(bool)
all_data = all_data[pre_filter_ix]
allX = allX[pre_filter_ix]
all_c = coord.SkyCoord(ra=all_data['ra']*u.degree, dec=all_data['dec']*u.degree)
In [ ]:
fig,axes = pl.subplots(1,2,figsize=(10,5),sharey=True)
for i in range(2):
ax = axes[i]
ax.plot(allX[:,i+1], allX[:,0], color='k', marker=',', ls='none', alpha=0.1)
for idx in np.argsort(all_data['non_ll'])[:100]:
thisX = allX[idx]
ax.plot(thisX[i+1], thisX[0], color='r', marker='o', ls='none')
axes[0].set_xlim(-0.5, 1.75)
axes[1].set_xlim(-0.5, 1.75)
ax.set_ylim(21, 14)
In [ ]:
iso_ll = all_data['iso_ll']
# iso_ll = all_data['xd_ll']
ll = (iso_ll - all_data['non_ll'])
weights = np.exp(ll)
# weights = np.exp(iso_ll-iso_ll.max())[filter_ix]
# filter_ix = non_ll > threshold
# filter_ix = np.ones_like(all_data['non_ll']).astype(bool)
# threshold total ll
threshold = -2
weights[ll > threshold] = np.exp(threshold)
# threshold control ll
# threshold = 10
# weights[all_data['non_ll'] < threshold] = np.exp(iso_ll[all_data['non_ll'] < threshold] - threshold)
In [ ]:
bins = np.linspace(-25, 25, 64)
pl.hist(iso_ll, bins=bins, alpha=0.5)
pl.hist(all_data['non_ll'], bins=bins, alpha=0.5)
pl.hist(ll, bins=bins, alpha=0.5)
pl.hist(np.log(weights), bins=bins, alpha=0.2)
# pl.axvline(threshold)
pl.yscale('log')
Plot CMD binned, weighted by weights
In [ ]:
cmd_bin_size = 0.02
g_bins = np.arange(allX[:,0].min(), allX[:,0].max(), step=cmd_bin_size)
gi_bins = np.arange(allX[:,1].min(), allX[:,1].max(), step=cmd_bin_size)
cmd_H,g_edges,gi_edges = np.histogram2d(allX[:,1], allX[:,0],
bins=(gi_bins, g_bins), weights=weights)
cmd_H_iso,g_edges,gi_edges = np.histogram2d(allX[:,1], allX[:,0],
bins=(gi_bins, g_bins), weights=np.exp(iso_ll-iso_ll.max()))
bg_H_iso,g_edges,gi_edges = np.histogram2d(allX[:,1], allX[:,0],
bins=(gi_bins, g_bins),
weights=np.exp(all_data['non_ll']-all_data['non_ll'].max()))
g_mesh,gi_mesh = np.meshgrid((g_bins[1:]+g_bins[:-1])/2, (gi_bins[1:]+gi_bins[:-1])/2)
In [ ]:
fig,axes = pl.subplots(1,3,figsize=(10,6),sharex=True,sharey=True)
ax = axes[0]
ax.pcolormesh(gi_mesh, g_mesh, cmd_H, cmap='Blues')
ax.plot(isoX[:,1], isoX[:,0]+fiducial_DM, marker=None)
ax.set_xlim(0, 0.75)
ax.set_ylim(21, 14)
ax.set_xlabel('$g-i$')
ax.set_ylabel('$i$')
ax = axes[1]
ax.pcolormesh(gi_mesh, g_mesh, cmd_H_iso, cmap='Blues')
ax.plot(isoX[:,1], isoX[:,0]+fiducial_DM, marker=None)
ax.set_xlabel('$g-i$')
ax = axes[2]
ax.pcolormesh(gi_mesh, g_mesh, bg_H_iso, cmap='Blues')
ax.plot(isoX[:,1], isoX[:,0]+fiducial_DM, marker=None)
ax.set_xlabel('$g-i$')
In [ ]:
sky_binsize = (6*u.arcmin).to(u.degree).value
sky_smooth = (6*u.arcmin).to(u.degree).value / sky_binsize
search_ra = ra
search_dec = dec
ra_bins = np.arange(search_ra.min(), search_ra.max()+sky_binsize, sky_binsize)
dec_bins = np.arange(search_dec.min(), search_dec.max()+sky_binsize, sky_binsize)
print("{} RA bins, {} Dec bins".format(len(ra_bins), len(dec_bins)))
In [ ]:
weights.shape
In [ ]:
# no_cluster_ix = all_c.separation(cluster_c) > 5*r_c
no_cluster_ix = np.ones(len(all_c)).astype(bool)
_data = all_data[no_cluster_ix]
H,ra_edges,dec_edges = np.histogram2d(_data['ra'], _data['dec'],
bins=(ra_bins,dec_bins), weights=weights[no_cluster_ix])
unw_H,ra_edges,dec_edges = np.histogram2d(_data['ra'], _data['dec'],
bins=(ra_bins,dec_bins))
ra_mesh,dec_mesh = np.meshgrid((ra_edges[1:]+ra_edges[:-1])/2, (dec_edges[1:]+dec_edges[:-1])/2)
H = H.T
unw_H = unw_H.T
# use a bivariate spline to smooth and divide out the background
spl = interpolate.SmoothBivariateSpline(ra_mesh.ravel(), dec_mesh.ravel(),
unw_H.ravel()/unw_H.sum(), kx=5, ky=5)
spl_unw_H = spl.ev(ra_mesh.ravel(), dec_mesh.ravel())
spl_unw_H = spl_unw_H.reshape(ra_mesh.shape)
# norm_H = H
norm_H = H / spl_unw_H
# norm_H = H / gaussian_filter(unw_H, 2*sky_smooth)
In [ ]:
# H_operation = lambda x: np.log(x)
# H_operation = lambda x: np.sqrt(x)
H_operation = lambda x: x
# H_operation = lambda x: x**2
In [ ]:
tmp = H_operation(norm_H.ravel())
bins = np.linspace(*scoreatpercentile(tmp, [1,99]), num=32)
pl.hist(tmp, bins=bins);
vmin,vmax = scoreatpercentile(tmp, [5,85])
pl.axvline(vmin, color='r')
pl.axvline(vmax, color='r')
In [ ]:
fig,axes = pl.subplots(3,1,figsize=(8,16),sharex=True,sharey=True)
ax = axes[0]
ax.pcolormesh(ra_mesh, dec_mesh, H_operation(norm_H),
cmap='Greys', vmin=vmin, vmax=vmax)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_c.to(u.degree).value,
edgecolor='r', facecolor='none')
ax.add_patch(pa)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_t.to(u.degree).value,
edgecolor='g', facecolor='none')
ax.add_patch(pa)
# ax.set_xlim(ra_mesh.max(), ra_mesh.min())
ax.set_xlim(ra_mesh.max(), 220.)
ax.set_ylim(dec_mesh.min(), dec_mesh.max())
ax.set_ylabel('Dec [deg]')
# ----
ax = axes[1]
ax.pcolormesh(ra_mesh, dec_mesh, gaussian_filter(H_operation(norm_H), sky_smooth),
cmap='Greys', vmin=vmin, vmax=vmax)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_c.to(u.degree).value,
edgecolor='r', facecolor='none')
ax.add_patch(pa)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_t.to(u.degree).value,
edgecolor='g', facecolor='none')
ax.add_patch(pa)
ax.set_ylabel('Dec [deg]')
# ----
ax = axes[2]
ax.pcolormesh(ra_mesh, dec_mesh, spl_unw_H,
cmap='Greys')#, vmin=vmin, vmax=vmax)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_c.to(u.degree).value,
edgecolor='r', facecolor='none')
ax.add_patch(pa)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_t.to(u.degree).value,
edgecolor='g', facecolor='none')
ax.add_patch(pa)
# ax.set_xlim(ra_mesh.max(), ra_mesh.min())
ax.set_xlim(ra_mesh.max(), 220.)
ax.set_ylim(dec_mesh.min(), dec_mesh.max())
ax.set_xlabel('RA [deg]')
ax.set_ylabel('Dec [deg]')
ax.set_aspect('equal')
fig.tight_layout()
In [ ]:
zoom_buffer = 2
fig,axes = pl.subplots(3,1,figsize=(6,14),sharex=True,sharey=True)
ax = axes[0]
ax.pcolormesh(ra_mesh, dec_mesh, H_operation(norm_H),
cmap='Greys', vmin=vmin, vmax=vmax)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_c.to(u.degree).value,
edgecolor='r', facecolor='none')
ax.add_patch(pa)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_t.to(u.degree).value,
edgecolor='g', facecolor='none')
ax.add_patch(pa)
ax.set_xlim(cluster_c.ra.degree+zoom_buffer, cluster_c.ra.degree-zoom_buffer)
ax.set_ylim(cluster_c.dec.degree-zoom_buffer, cluster_c.dec.degree+zoom_buffer)
# ax.set_xlabel('RA [deg]')
ax.set_ylabel('Dec [deg]')
ax.set_aspect('equal')
# ----
ax = axes[1]
ax.pcolormesh(ra_mesh, dec_mesh, gaussian_filter(H_operation(norm_H), sky_smooth),
cmap='Greys', vmin=vmin, vmax=vmax)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_c.to(u.degree).value,
edgecolor='r', facecolor='none')
ax.add_patch(pa)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_t.to(u.degree).value,
edgecolor='g', facecolor='none')
ax.add_patch(pa)
ax.set_ylabel('Dec [deg]')
# ----
ax = axes[2]
ax.pcolormesh(ra_mesh, dec_mesh, spl_unw_H,
cmap='Greys')#, vmin=vmin, vmax=vmax)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_c.to(u.degree).value,
edgecolor='r', facecolor='none')
ax.add_patch(pa)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_t.to(u.degree).value,
edgecolor='g', facecolor='none')
ax.add_patch(pa)
ax.set_xlim(cluster_c.ra.degree+zoom_buffer, cluster_c.ra.degree-zoom_buffer)
ax.set_ylim(cluster_c.dec.degree-zoom_buffer, cluster_c.dec.degree+zoom_buffer)
ax.set_xlabel('RA [deg]')
ax.set_ylabel('Dec [deg]')
ax.set_aspect('equal')
fig.tight_layout()
In [ ]:
sm_H = gaussian_filter(H_operation(norm_H), sky_smooth)
levels = scoreatpercentile(norm_H.ravel(), [55,60,70,80,90,95,96,97,98,99])
# levels = scoreatpercentile(norm_H.ravel(), np.linspace(85,100,16))
fig,ax = pl.subplots(1,1,figsize=(6,6),sharex=True,sharey=True)
ax.contour(ra_mesh, dec_mesh, sm_H, levels=levels, colors='k')
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_c.to(u.degree).value,
edgecolor='r', facecolor='none')
ax.add_patch(pa)
pa = mpl.patches.Circle((cluster_c.ra.degree, cluster_c.dec.degree), radius=r_t.to(u.degree).value,
edgecolor='g', facecolor='none')
ax.add_patch(pa)
ax.set_xlim(cluster_c.ra.degree + zoom_buffer, cluster_c.ra.degree - zoom_buffer)
ax.set_ylim(cluster_c.dec.degree - zoom_buffer, cluster_c.dec.degree + zoom_buffer)
ax.set_xlabel('RA [deg]')
ax.set_ylabel('Dec [deg]')
ax.set_aspect('equal')
In [ ]:
In [ ]:
In [ ]:
In [ ]:
def rot_matrix(t, v):
m1 = 1 - np.cos(t)
R = np.array([[np.cos(t) + v[0]**2*m1, v[0]*v[1]*m1 - v[2]*np.sin(t), v[1]*np.sin(t)+v[0]*v[2]*m1],
[v[2]*np.sin(t)+v[0]*v[1]*m1, np.cos(t)+v[1]**2*m1, -v[0]*np.sin(t)+v[1]*v[2]*m1],
[-v[1]*np.sin(t)+v[0]*v[2]*m1, v[0]*np.sin(t)+v[1]*v[2]*m1, np.cos(t)+v[2]**2*m1]])
return R
Fit a polynomial surface
In [ ]:
# vec = cluster_c.cartesian.xyz.value
# c = coord.ICRS(ra=all_data['ra']*u.degree, dec=all_data['dec']*u.degree)
# def surface_model(xdata, *p):
# theta,a,b,c1,c2,c3,c4,c5 = p
# R = rot_matrix(theta, vec)
# rot_xyz = np.array([R.dot(w).tolist() for w in c.cartesian.xyz.value.T])
# sph = coord.CartesianRepresentation(rot_xyz.T*u.one)\
# .represent_as(coord.UnitSphericalRepresentation)
# x = sph.lon.degree
# return a*np.exp(-b*x) * polyval(x, c=[c1, c2, c3, c4, c5])
In [ ]:
# curve_fit(surface_model, )
Transform to "stream coordinates"
In [ ]:
import astropy.coordinates as coord
In [ ]:
vv = coord.ICRS(ra=230*u.degree, dec=-21*u.degree).cartesian.xyz.value
R = rot_matrix((52+180)*u.degree, vv)
c = coord.ICRS(ra=all_data['ra']*u.degree, dec=_dec*u.degree)
new_xyz = np.array([R.dot(w).tolist() for w in c.cartesian.xyz.value.T])
In [ ]:
sph = coord.CartesianRepresentation(new_xyz.T*u.one)\
.represent_as(coord.UnitSphericalRepresentation)
In [ ]:
pl.plot(sph.lon.degree, sph.lat.degree, ls='none', marker=',', alpha=0.051)
In [ ]:
phi1, phi2 = sph.lon.degree, sph.lat.degree
sph_H,phi1_edges,phi12_edges = np.histogram2d(phi1, phi2,
bins=(ra_bins+5,dec_bins),
weights=weights, normed=True)
phi1_mesh,phi2_mesh = np.meshgrid((phi1_edges[1:]+phi1_edges[:-1])/2,
(phi12_edges[1:]+phi12_edges[:-1])/2)
sph_H = sph_H.T
In [ ]:
pl.figure(figsize=(8,8))
pl.pcolormesh(sph_H, cmap='Greys', vmin=vmin, vmax=vmax)
x_lims = (70, 145)
for l in x_lims:
pl.axvline(l, color='r')
In [ ]:
sph_H_cut = sph_H[:, slice(*x_lims)]
In [ ]:
pl.figure(figsize=(8,8))
pl.pcolormesh(sph_H_cut, cmap='Greys', vmin=vmin, vmax=vmax)
pl.gca().set_aspect('equal')
slc_lo = 34
slc_hi = 44
pl.axhline(slc_lo, color='r')
pl.axhline(slc_hi, color='r')
In [ ]:
from scipy.optimize import curve_fit
from scipy.stats import scoreatpercentile
from numpy.polynomial.polynomial import polyval
In [ ]:
fit_y = np.mean(sph_H_cut[slc_lo:slc_hi], axis=0)
fit_x = np.arange(fit_y.size)
y_err = 0.05 * fit_y.mean() # 5% of mean
# just exponential
# def model(x, a, b):
# return a*np.exp(-b*x)
# p0 = (scoreatpercentile(fit_y, 95), 0.78E-2) # 0.78 - BY EYE fit
# exponential * polynomial
def model(x, a, b, c1, c2, c3, c4):
return a*np.exp(-b*x) * polyval(x, c=[c1, c2, c3, c4])
p0 = (scoreatpercentile(fit_y, 95), 1E-2, 0, 0, 0, 0)
In [ ]:
p_opt, pcov = curve_fit(model, xdata=fit_x, ydata=fit_y, p0=p0, sigma=y_err)
In [ ]:
pl.figure(figsize=(10,8))
for l,r in zip(range(0,sph_H_cut.shape[0],10), range(10,sph_H_cut.shape[0]+10,10)):
pl.plot(np.mean(sph_H_cut[l:r], axis=0), marker=None)
_x = np.linspace(0,sph_H_cut.shape[1],256)
pl.plot(_x, model(_x, *p_opt), marker=None, lw=3, ls='--')
# pl.plot(_x, p0[0]*np.exp(-5E-3*(_x)), marker=None, lw=3, ls='--')
# pl.plot(_x, model(_x, 0.01481421, 9.85E-3, 0.1), marker=None, lw=3, ls='--')
pl.ylim(8E-4,8E-2)
# pl.xlim(1,100)
pl.yscale('log')
# pl.xscale('log')
In [ ]:
_gridx = np.arange(sph_H_cut.shape[0])
_gridy = np.arange(sph_H_cut.shape[1])
_xx,_ = np.meshgrid(_gridy,_gridx)
model_im = model(_xx, *p_opt)
In [ ]:
div_rot_H = (sph_H_cut - model_im)
In [ ]:
pl.figure(figsize=(10,8))
for l,r in zip(range(0,sph_H_cut.shape[0],10), range(10,sph_H_cut.shape[0]+10,10)):
pl.plot(np.mean(div_rot_H[l:r], axis=0), marker=None)
pl.ylim(-0.01,0.01)
# pl.xlim(1,100)
# pl.yscale('log')
# pl.xscale('log')
In [ ]:
drH_rav = div_rot_H.ravel()
pl.hist(drH_rav, bins=np.linspace(-0.01,0.05,64));
dr_vmin,dr_vmax = scoreatpercentile(drH_rav, [15,85])
pl.axvline(dr_vmin, color='r')
pl.axvline(dr_vmax, color='r')
In [ ]:
pl.figure(figsize=(7,7))
pl.pcolormesh(np.flipud(div_rot_H), cmap='Greys', vmin=dr_vmin, vmax=0)
pl.gca().set_aspect('equal')
In [ ]:
pl.figure(figsize=(7,7))
H_sigma = (6*u.arcmin / bin_size).decompose().value
pl.pcolormesh(gaussian_filter(np.flipud(div_rot_H), H_sigma),
cmap='Greys', vmin=dr_vmin, vmax=dr_vmax)
pl.gca().set_aspect('equal')
In [ ]:
sub = (slice(60,100), slice(30,120))
In [ ]:
pl.figure(figsize=(7,7))
H_sigma = (6*u.arcmin / bin_size).decompose().value
pl.pcolormesh(gaussian_filter(np.flipud(div_rot_H)[sub], H_sigma),
cmap='Greys', vmin=dr_vmin, vmax=dr_vmax)
pl.gca().set_aspect('equal')
In [ ]: